%% fig3a: Lyapunov 指数谱（Benettin 方法，parfor 并行计算）
clc; clear; close all;

%% —— 参数设置 —— 
g11        = 1.4;  
g12        = -2;
g21        = 1.4;  
g22        = 1.6;

x0         = [0.01; 0; 0];    % 初始状态
dt         = 0.01;           
t_end      = 300;             % 总积分时间
renormStep = 50;               % 每隔 renormStep 步正交化一次

k_list     = linspace(0,3,300);
LEs_all    = zeros(3, numel(k_list));

%% —— 并行计算 Lyapunov 指数谱 —— 
parfor idx = 1:numel(k_list)
    k = k_list(idx);

    % 系统方程句柄
    f = @(x) MBNHNN(x, k, g11, g12, g21, g22);
    % 解析雅可比句柄
    J = @(x) Jaco(x, k, g11, g12, g21, g22);

    % Benettin 方法计算全谱
    LEs_all(:,idx) = Benettin(f, J, x0, dt, t_end, renormStep);
end

%% —— 绘图 —— 
figure('Color','w');
plot(k_list, LEs_all(1,:), 'r-', 'LineWidth',1.5); hold on;
plot(k_list, LEs_all(2,:), 'g-', 'LineWidth',1.5);
plot(k_list, LEs_all(3,:), 'b-', 'LineWidth',1.5);
xlabel('k', 'FontSize',14);
ylabel('LEs', 'FontSize',14);
legend('\lambda_1','\lambda_2','\lambda_3','Location','Best');
grid on;
ylim([-0.7, 0.3]);


